home *** CD-ROM | disk | FTP | other *** search
/ CD ROM Paradise Collection 4 / CD ROM Paradise Collection 4 1995 Nov.iso / science / fastmap.zip / HOMOGC.C < prev    next >
C/C++ Source or Header  |  1991-12-09  |  7KB  |  245 lines
  1. /***************************************
  2.  *                                     *
  3.  *             HOMOG                   *
  4.  *                                     * 
  5.  *    Copyright Dr. Jurg Ott 1985      *
  6.  *                                     *
  7.  *    transcribed from Fortran to C    *
  8.  *    by Dr. David Curtis 1988         *
  9.  ***************************************/
  10.  
  11. /*
  12. In:
  13. Jurg Ott 1985 Analysis of Human Genetic Linkage, John Hopkins University Press 
  14.  
  15. To carry out the A-test (Ott,Ann Hum Genet 47:311-320,1983) estimating the 
  16. proportion alpha of "linked" families jointly with the recombination fraction 
  17. theta in these families. Input (free format) and output are interactive or on 
  18. files.
  19.  
  20.                       Input information
  21.  
  22. NR = number of lod scores per family (do not include z=0 at rec. freq. 0.5)
  23. ITAB = 1 if lod table is desired, = 0 for test and estimates only
  24. NR rec. freq. (theta) values in increasing order, e.g. .01 .05  ...
  25. Then input the following information for each family:
  26. ISW = 2 if number of recombinants and nonrecombinants are given for next family
  27.     = 1 if lod scores follow for next family
  28.     = 0 to finish entering family information
  29. V[i] or RK1 RK2 = lod scores or number of recombinants for family. Note that 
  30. lods smaller than -30, e.g. -35, are taken to represent minus infinity and 
  31. will appear on output as -99. 
  32. */
  33.  
  34. #include <stdio.h>
  35. #include <math.h>
  36. FILE *fi,*fo;
  37. #define MAXL 20              /* maximum number of lod scores */
  38. #define CONS .30103
  39. #define DLIM -30.
  40. #define ELIM 1.0E-30
  41. #define MAXF 50              /* maximum number of families   */
  42. #define STEP .05
  43. int ISW,ITAB,INF[MAXL],NF,i,j,l,NR,NR1;
  44. double RL[MAXF][MAXL],V[MAXL],VS[MAXL+1],THET[MAXL+1];
  45. double AL,R;
  46.  
  47. double pnorm();
  48. void wr_VS(),comp_res(),get_data();
  49. FILE *get_file();
  50.  
  51. double pnorm(x2)
  52. double x2;
  53.  {
  54.  double t;
  55.  if (x2<=0) return .5;
  56.  t=1/(1+.33267*sqrt(x2));
  57.  return .39894228*exp(-.5*x2)*t*(.4361836+t*(.937298*t-.1201676));
  58.  }
  59.  
  60. void wr_VS()
  61.   {
  62.   int l;
  63.   fprintf(fo," %5.2lf  ",AL);
  64.   for (l=1;l<=NR1;++l) 
  65.        {
  66.        fprintf(fo,"%10.4lf",VS[l]);
  67.        if (l==10) fprintf(fo,"\n        ");
  68.        }
  69.   fprintf(fo,"\n");
  70.   }
  71.  
  72. FILE *get_file(mode)
  73. char *mode;
  74.  {
  75.  char line[240],s[20];
  76.  FILE *fp;
  77.  do {
  78.  *s=*line='\0';
  79.  printf("\nEnter %sput file name (just [ENTER] for keyboard) ---> ",mode[0]=='r'?"in":"out") ;
  80.  gets(line);
  81.  if (sscanf(line," %s",s)!=1) return mode[0]=='r'?stdin:stdout;
  82.     } while ((fp=fopen(s,mode))==0 && (printf("\nCould not open %s\n",s),1));
  83.  return fp;
  84.  }
  85.  
  86. main()
  87.  {
  88.  fi=get_file("r");
  89.  fo=get_file("w");
  90.  fprintf(fo,"\n  HOMOGENEITY TESTS\n");
  91.  for (;;)
  92.  {
  93.  get_data();
  94.  comp_res();
  95.  }
  96.  return 0;
  97.  }
  98.  
  99. void get_data()
  100.  {
  101.  double RKK,RK1,RK2;
  102.  printf("\n Enter number of lod scores per family (0 to exit program) ---> ");
  103.  fscanf(fi," %d",&NR);
  104.  if (!NR) exit (0);
  105.  if (NR>MAXL) {fprintf(stderr,"Number of lod scores too large\n");exit(1);}
  106.  NR1=NR+1;
  107.  printf("\nLod tables desired? (1=Yes, 0=No) ---> ");
  108.  fscanf(fi," %d",&ITAB);
  109.  printf("\nEnter %d rec. fr. values:\n",NR);
  110.  for (i=1;i<=NR;++i)
  111.     {
  112.     fscanf(fi," %lf",&THET[i]);
  113. #if 0
  114.     if (THET[i]>1) 
  115.         {printf(stderr,"Error: Rec. fr. values must be less than 1\n");exit(3);}
  116.     /* Allow user to enter map distances or whatever instead of thetas
  117.        - A-test should be just as valid, provided lod scores are entered
  118.        and not number of recombinant/non-recombinants.
  119.     */
  120. #endif
  121.     INF[i]=VS[i]=0;
  122.     }
  123.  THET[NR1]=0.5;
  124.  VS[NR1]=0;
  125.  for (NF=0;;)
  126.   {
  127.   puts("\nNew family? (2=yes, with no. of rec. and nonrec.;");
  128.   printf(" 1=yes, with lod scores; 0=no, carry out test) ---> ");
  129.   fscanf(fi," %d",&ISW);
  130.   if (!ISW) return;
  131.   ++NF;
  132.   if (NF>MAXF) {fprintf(stderr,"\nNumber of families too large!\n");exit(2);}
  133.   if (ISW!=1)
  134.     {
  135.     printf("\nEnter number of recombinants and number of non-rec. ---> ");
  136.     fscanf(fi," %lf %lf",&RK1,&RK2);
  137.     RKK=CONS*(RK1+RK2);
  138.     for (i=1;i<=NR;++i)
  139.       {
  140.       R=THET[i];
  141.       V[i]=(RK1&&!R)?-90:RKK+RK1*log10(R)+RK2*log10(1-R);
  142.       }
  143.     }
  144.   else
  145.     {
  146.     printf("\nEnter %d lod scores  --->  ",NR);
  147.     for (i=1;i<=NR;++i) 
  148.       fscanf(fi,"%lf",&V[i]);
  149.     }
  150.   for (i=1;i<=NR;++i)
  151.     {
  152.     VS[i]+=(R=V[i]);
  153.     if (R<DLIM)
  154.        {RL[NF][i]=0;INF[i]=1;}
  155.     else RL[NF][i]=pow(10.,R);
  156.     }
  157.   }
  158.  }
  159.  
  160. void comp_res()
  161.  {
  162.  int JM=NR1,JT;
  163.  double VM=0,VT,SUM;
  164.  double ALT,CC1,C2,CCC2,E1,E2,E3,AL1;
  165.  for (i=1;i<=NR;++i)
  166.   {
  167.   if (INF[i]) VS[i]=-99;
  168.   else if (VS[i]>VM)
  169.      {
  170.      VM=VS[i];
  171.      JM=i;
  172.      }
  173.   }
  174.   CC1=4.60517*VM;
  175.   E2=pnorm(CC1);
  176.   JT=JM;
  177.   ALT=1;
  178.   VT=VM;
  179.   AL=1;
  180.   if (ITAB)
  181.     {
  182.     fprintf(fo,"\n Lod table\n\n Alpha\n");
  183.     wr_VS();
  184.     }
  185.   while ((AL-=STEP)>.01)
  186.     {
  187.     AL1=1-AL;
  188.     for (l=1;l<=NR;++l)
  189.       {
  190.       SUM=1;
  191.       for (i=1;i<=NF;++i) SUM*=AL*RL[i][l]+AL1;
  192.                               /* Note: in original Fortran listing */
  193.                               /* this line should read:            */
  194.                               /*    SUM=SUM*(AL*RL(I,L)+AL1)       */
  195.                               /* and not as printed                */
  196.       if (SUM<ELIM) SUM=-99;
  197.       else 
  198.         {
  199.         SUM=log10(SUM);
  200.         if (SUM>VT)
  201.          {
  202.          VT=SUM;
  203.          ALT=AL;
  204.          JT=l;
  205.          }
  206.         }
  207.       VS[l]=SUM;
  208.       }
  209.     if (ITAB) wr_VS();
  210.     }
  211.  if (ITAB)
  212.    {
  213.    AL=0;
  214.    for(l=1;l<=NR;++l) VS[l]=0;
  215.    wr_VS();
  216.    fprintf(fo,"\n        ");for(j=1;j<=NR1;++j) fprintf(fo,"%10.4lf",THET[j]);
  217.    fprintf(fo,"\n                                        Theta\n");
  218.    }
  219.  if (JT==NR1) C2=CCC2=ALT=0;
  220.  else
  221.    {
  222.    C2=4.60517*(VT-VM);
  223.    CCC2=4.60517*VT;
  224.    }
  225.  E1=pnorm(C2);
  226.  E3=0.5*exp(-.5*CCC2);
  227.  fprintf(fo,"\n Hypotheses  Max. lod  Alpha  Theta\n");
  228.  fprintf(fo," H2%18.3lf%7.2lf%7.2lf\n",VT,ALT,THET[JT]);
  229.  fprintf(fo," H1%18.3lf    (1)%7.2lf\n H0              (0)     (0)   (0.5)",VM,THET[JM]);
  230.  fprintf(fo,"\n Components of chi-square\n\n Source                  ");
  231.  fprintf(fo,"d.f.   chi-square  p-value\n");
  232.  fprintf(fo," H2 vs. H1  Heterogeneity  1%13.3lf%9.4lf\n",C2,E1);
  233.  fprintf(fo," H1 vs. H0  Linkage        1%13.3lf%9.4lf\n",CC1,E2);
  234.  fprintf(fo," H2 vs. H0  Total          2%13.3lf%9.4lf\n",CCC2,E3);
  235.  fprintf(fo,"\n\n Family nr.  Posterior prob. of linkage\n");
  236.  AL=1-ALT;
  237.  for (i=1;i<=NF;++i)
  238.    {
  239.    R=ALT*RL[i][JT];
  240.    R/=(R+AL);
  241.    fprintf(fo," %8d%22.4lf\n",i,R);
  242.    }
  243.  }
  244.  
  245.